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ABSTRACT 

The  detection  performance  of  maritime  radars  is  usually  limited  by  sea  clutter.  The  K  distribution 
is  a  well  established  statistical  model  of  sea  clutter  which  is  widely  used  in  performance 
calculations.  There  is  no  closed  form  solution  for  the  probability  of  detection  in  K-distributed 
clutter,  so  numerical  methods  are  required.  The  K  distribution  is  a  compound  model  which 
consists  of  Gaussian  speckle  modulated  by  a  slowly  varying  mean  level,  this  local  mean  being 
gamma  distributed.  A  series  solution  for  the  probability  of  detection  in  Gaussian  noise  is 
integrated  over  the  gamma  distribution  for  the  local  clutter  power.  Gauss-Laguerre  quadrature  is 
used  for  the  integration,  with  the  nodes  and  weights  calculated  using  matrix  methods,  so  that  a 
general  purpose  numerical  integration  routine  is  not  required.  The  method  is  implemented  in 
Matlab  and  compared  with  an  approximate  solution  based  on  lookup  tables.  The  solution 
described  here  is  slower,  but  more  accurate  and  more  flexible  in  that  it  allows  for  a  wider  range  of 
target  fluctuation  models. 
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Calculation  of  Radar  Probability  of  Detection  in 
K-Distributed  Sea  Clutter  and  Noise 


Executive  Summary 

Modelling  of  the  radar  returns  from  the  sea  is  required  for  operations  analysis  of  maritime 
patrol,  in  order  to  calculate  probabilities  of  detection  for  targets  of  interest.  The  detection 
performance  of  maritime  radars  is  usually  limited  by  sea  clutter.  A  statistical  model  of  the 
clutter  is  normally  used,  and  the  K  distribution  has  become  standard  for  this  purpose. 
There  is  no  closed  form  solution  for  the  probability  of  detection  in  K-distributed  clutter,  so 
numerical  methods  are  required.  A  method  for  calculation  of  the  probability  of  detection 
in  K-distributed  clutter  and  noise  is  described,  and  compared  with  other  approaches.  It  is 
faster  than  direct  numerical  integration,  and  more  flexible  and  more  accurate,  but  slower, 
than  an  approximate  method  based  on  interpolation.  Faster  evaluation  is  highly  desirable 
for  simulation  models  used  in  operations  research.  The  method  described  here  could  be 
used  to  create  tables  for  interpolation,  which  should  produce  the  appropriate  mix  of  speed 
and  accuracy  for  operations  research  simulation  models. 
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1.  Introduction 


Modelling  of  the  radar  returns  from  the  sea  is  required  for  operations  analysis  of  maritime 
patrol,  in  order  to  calculate  probabilities  of  detection  for  targets  of  interest.  In  this  context  the 
backscatter  from  the  sea  is  referred  to  as  sea  clutter.  A  statistical  model  of  the  clutter  is 
normally  used,  and  the  K  distribution  has  become  standard  for  this  purpose  [1],  Sea  clutter 
consists  of  a  rapidly  varying  speckle  component  and  an  underlying  mean  amplitude  which 
varies  more  slowly.  The  K  distribution  provides  a  compound  representation  which  includes 
both  components.  Thermal  noise,  which  has  Gaussian  statistics,  must  also  be  considered  in 
calculating  probabilities  of  detection.  Fluctuations  in  the  target  return  also  need  to  be  taken 
into  account,  and  this  is  normally  done  with  statistical  models  based  on  the  gamma 
distribution. 

The  calculation  of  the  probability  of  detection  in  K-distributed  clutter  and  noise  is  quite 
difficult.  An  approximate  method  was  developed  by  Watts  and  Wicks  [2,3].  This  method  is 
fast,  because  the  probability  of  detection  is  obtained  by  linear  interpolation  in  a  table, 
following  some  preliminary  calculations  which  need  only  be  done  once  for  each  set  of 
parameters.  Tlowever,  the  result  is  only  approximate,  and  the  necessary  tables  of  coefficients 
[3]  are  only  available  for  the  Swerling  1  and  2  target  fluctuation  models,  as  well  as  a  non¬ 
fluctuating  target  (Swerling  0).  Alternatively,  the  probability  of  detection  can  be  calculated  by 
numerical  integration  over  the  K  distribution.  This  will  give  accurate  results,  but  it  is  slow, 
and  therefore  impractical  for  use  in  simulations  where  the  probability  of  detection  must  be 
updated  frequently.  Another  approach  is  described  in  [1],  and  further  explored  in  this  report, 
in  which  the  K  distribution  is  separated  into  its  components  and  only  the  integration  over  the 
local  clutter  power  is  carried  out  numerically. 


2.  Probability  of  Detection 


2.1  The  K  Distribution 

The  K  distribution  for  the  clutter  intensity  can  be  expressed  as 

/•  00 

P(z)  =  J  P(z  x)Pc(x)dx 

(1) 

where 

P(z  x)  =  —  exp(-z/x) 

X 

(2) 

and 

hv  xv' 1 

pc  0)=  exp(  bx) 

(3) 

T(v) 

Equation  (3)  is  the  Gamma  distribution  for  the  local  clutter  power  x.  The  mean  power  is 
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pc=(x)  =  v/b 

Equation  (1)  evaluates  to 

p{z)  =  ^-(4b~z)V'X  Kv_\ 2Vfe) 


(4) 

(5) 


The  modified  Bessel  function  K  in  equation  (5)  gives  the  distribution  its  name.  However,  for 
calculating  the  probability  of  detection  it  is  advantageous  to  keep  the  speckle  and  modulation 
components  separate,  assuming  that  the  local  clutter  power  x  remains  constant  over  the  beam 
dwell  time. 


2.2  Probability  of  Detection  Calculation 

The  probability  of  detection  in  Gaussian  noise  and  speckle  is 

/•  oo 

Pd{Y\x,N)  =  \YPR{p\s,N)d/u  (6) 


The  threshold  Y  is  normalised  by  the  noise  and  local  clutter  power  x. 

The  sum  of  N  radar  returns,  normalised  by  the  noise  and  local  clutter  power  is 


1  N 

X  +  Pn  /= 1 


(7) 


The  sum  of  the  target  powers  from  the  N  pulses  is 


1  N 

X+Pn  i= 1 


(8) 


The  probability  density  function  for  the  sum  ji  of  N  radar  returns  from  a  fixed  target  in 
Gaussian  noise  is  the  multilook  Rice  distribution. 


Pr(p\s,N)  =  \^ 
\s  J 


\(N- 1)/2 


-(/j+s)  i 


-i(2V/^) 


-u  N+i-l  -s  i 

e  p  e  s 


rj(N  +  i -1)1  /! 


(9) 


z 


The  series  form  of  equation  (9)  is  obtained  using  the  series  expansion  for  the  Bessel  function 
[4,  §8.445]: 


'„(*)  =  £  T 


1 


ttiT(n  +  i  +  \) 


(  z\n+2i 


(10) 
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The  integral  in  equation  (6)  can  be  evaluated  to  give  the  double  sum  formula  for  the 
probability  of  detection  for  a  fixed  target  [5,6,7]: 


PR{M\s,N)dju  =  YJ 


00  e~ssl 
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,=o  /!  {N  +  i-\)Vr 

I" 


f  00  -//  N+i-l  i 

J  y  e  A  d/u 


°°  a-s  J  N+i-l  yj 

—  &  b  _y  ^  I 
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YJ 
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The  manipulation  of  the  series  in  equation  (11)  makes  use  of  the  result 

®  p~sY 


i= 0 


Z! 


(ii) 


(12) 


Target  fluctuations  are  modelled  using  the  Gamma  distribution  for  the  sum  of  target  returns  s: 


PT(s  |  S',  At)  = 


T(k) 


k 

J 


-ks/S 


where 


n(a2) 

x+p„ 


(13) 


(14) 


The  Gamma  distribution  encompasses  all  the  Marcum-Swerling  and  Weinstock  target  models 
[8], [9,  Section  B]  as  shown  in  the  table  below: 


Weinstock 

0  <  k  <  1 

Swerling  1 

k  =  1 

Swerling  2 

k  =  N 

Swerling  3 

k  =  2 

Swerling  4 

Z 

(N 

II 

Swerling  0 

k  — »  oo 

A  non-fluctuating  target  (Swerling  0)  corresponds  to  k  — *•  oo  or  s  =  S. 
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The  probability  of  detection  for  a  fluctuating  target  is  [9,  Section  B] 

/•co  /•  co 

Pd(Y  |  x,S,k,N )  =  j  Pr(s  |  PR{ju  |  s,N)djU  ds 

YJ 


N- 1 

j=0 
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7!  h  j 
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(15) 


f  k  > 

k 

(  s  1 

A 

[s  +  k  j 

[s  +  kj 

J 

Numerical  calculation  of  the  probability  of  detection  from  equations  (11)  and  (15)  is  not 
straightforward,  due  to  underflow  and  loss  of  significance.  Appendix  A  gives  some  insight 
into  the  nature  of  the  series  in  these  equations,  and  shows  why  it  is  not  always  possible  to 
begin  the  summation  with  the  first  term,  even  with  high  precision  arithmetic.  Shnidman 
[5,6,7,9,10]  has  developed  methods  to  carry  out  these  calculations  efficiently,  and  coded  the 
resulting  algorithms  in  Matlab  L  Chernoff  bounds  (Appendix  B)  are  used  to  avoid 
unnecessary  calculations. 


2.3  Integration  over  the  Clutter  Power 

The  probability  of  detection  in  K-distributed  clutter  with  a  single  pulse  threshold  1 /„  and  local 
clutter  power  x  is 

PAyn^)  =  j0  Pd(Y\x,N)Pc(x)dx  (16) 


where  Pd(Y  \  x)  is  calculated  as  for  detection  in  noise,  equation  (6),  and  Pc(x)  is  the  Gamma 
distribution  (3).  The  threshold  yn  is  normalised  by  the  sum  of  the  noise  power  p>„  and  the  mean 
clutter  power  pc.  At  first  sight,  numerical  integration  over  the  local  clutter  power  seems 
unlikely  to  provide  a  fast  and  efficient  method  for  calculating  the  probability  of  detection. 
However,  if  we  make  the  substitution  t  =  bx  in  the  Gamma  distribution  (3)  this  becomes 


Pc(x)dx  =  Pc(t)dt  = 


Y(v) 


-dt 


(17) 


In  terms  of  the  mean  clutter  power  pCr  t  =  vx/ pc.  Equation  (16)  becomes 


PAy„N) 


/•CO  , 

- [  P,(Y\t,  N)r'e-'dt 

F(v)Jo 


(18) 


This  integral  is  readily  evaluated  using  Gauss-Laguerre  quadrature  [11,  Section  4.5],  The 
Laguerre  polynomials  are  generated  from  a  recurrence  relation,  and  the  nodes  and  weights 
are  calculated  from  the  eigenvalues  and  eigenvectors  of  a  symmetric  tridiagonal  matrix. 
Gautschi  [12]  has  written  Matlab  code  to  do  this.  In  most  cases  sufficient  accuracy  in  Pd  is 
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obtained  with  only  10  quadrature  points.  The  exception  is  low  thresholds  and  v  <  1,  where 
many  more  points  are  needed  to  cope  with  the  singularity  at  the  origin.  However,  this 
situation  is  not  of  much  practical  importance,  because  a  high  threshold  is  needed  to  achieve  a 
satisfactory  probability  of  false  alarm  when  v  is  small.  The  maximum  value  of  v  for  which  this 
method  of  integration  can  be  used  in  Matlab  is  171,  because  it  requires  T (v)  to  be  less  than  the 
largest  floating  point  number  which  can  be  represented,  i.e.  T (v)  <  10J°8.  As  v  — >  oo  the  clutter 
becomes  noise  like,  so  integration  over  the  Gamma  distribution  is  not  required.  There  is  a  very 
small,  but  perceptible,  difference  between  the  probability  of  detection  calculated  for  v  =  170 
and  the  Rayleigh  limit  where  the  clutter  is  treated  as  noise. 

2.4  Normalisation 

The  threshold  Y  is  normalised  by  the  noise  and  local  clutter  power.  In  terms  of  a  single  pulse, 
unnormalised  threshold  y 


y_  Ny 

X+Pn 


(19) 


We  would  like  to  express  this  in  terms  of  the  clutter  to  noise  ratio  CNR  and  the  integration 
variable  t. 


y  N{pc+pn) 

Pc+Pn  tP‘  +  pn 

V 

N(\  +  CNR ) 

y  n  f 

1  +  —  CNR 


(20) 


where  y«  is  the  single  pulse  threshold  normalised  by  the  mean  clutter  plus  noise  power. 


In  the  same  way  we  would  like  to  express  the  signal  S  in  terms  of  the  signal  to  noise  ratio 
SNR,  the  clutter  to  noise  ratio  CNR  and  the  integration  variable  t.  From  equation  (14), 


n(a2) 

X  +  Pn 


N  SNR 

1 +  -CNR 

v 


(21) 
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For  the  clutter  only  case  we  have  p„  =  0  or  CNR  — >  oo  so 


y  =  y, 


Nv 

t 


Nv 

and  5  = - SCR 

t 


(22) 


where  SCR  is  the  signal  to  clutter  ratio. 

2.5  Probability  of  False  Alarm 

The  probability  of  false  alarm  is  obtained  from  equation  (11)  with  s  =  0  (note  that  0°  =  1) 


AM  yj 

Pfa(Y  \x,N)  =  Y,e-Y  — 


r(N,Y) 

F(A0 


(23) 


In  terms  of  the  Matlab  incomplete  gamma  function,  Pf  =  1  -  gammainc(Y,N). 

In  order  to  calculate  the  probability  of  detection  for  a  given  probability  of  false  alarm,  the 
threshold  must  be  determined  from  the  probability  of  false  alarm.  This  can  be  done  using  the 
Matlab  root  finding  function  fzero.  This  function  requires  either  a  starting  estimate  for  the 
root,  or  an  interval  which  brackets  the  root.  In  the  first  case  fzero  searches  for  the  starting 
interval  itself,  which  tends  to  cause  it  to  fail,  because  the  probability  of  false  alarm  is 
undefined  for  negative  thresholds,  and  underflows  to  zero  for  large  positive  thresholds.  A 
starting  interval  of  0.1  to  1100  for  yn  was  found  to  be  satisfactory.  Integration  over  the  local 
clutter  power  does  not  impede  the  calculation  much,  because  the  nodes  and  weights  depend 
only  on  v,  so  they  don't  need  to  be  recalculated  during  the  root  finding  process. 


In  the  case  of  a  single  pulse  return  from  clutter,  the  integral  over  the  clutter  power  can  be 
evaluated  to  give 


P fa  (y n)  — 


1 


r(v) 

2(vy„r/2 

r(v) 


r  00  i 

J0  e-'dt 

(2CC,) 


K 


(24) 


This  solution  is  used  in  Section  3.4  to  assess  the  accuracy  of  the  numerical  integration  with 
different  numbers  of  quadrature  points. 
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3.  Sample  Calculations 

3.1  Probability  of  False  Alarm 

Figure  1  shows  the  probability  of  false  alarm  for  a  single  pulse  return  from  K-distributed 
clutter,  for  different  values  of  the  shape  parameter  v.  Figure  1  reproduces  Figure  8.6  from 
reference  [1], 


Figure  1  Probability  of  false  alarm  for  a  single  pulse  return  from  K-distributed  clutter,  for  different 

values  of  the  shape  parameter  v 

3.2  Probability  of  Detection 

Figure  2  shows  the  probability  of  detection  for  a  single  pulse  return  from  a  fixed  target  in 
K-distributed  clutter  with  v  =  10,  for  different  probabilities  of  false  alarm.  Figure  2  reproduces 
the  top  graph  in  Figure  8.11  from  reference  [1],  Figure  3  shows  the  probability  of  detection  for 
the  incoherent  sum  of  10  pulses  from  a  Swerling  1  target  in  K-distributed  clutter  with  v  =  0.1, 
for  different  probabilities  of  false  alarm.  Figure  3  reproduces  the  bottom  graph  in  Figure  8.14 
from  reference  [1], 
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Figure  2  Probability  of  detection  for  a  single  pulse  return  from  a  fixed  target  in  K-distributed  clutter 
with  v  =  10,  for  different  probabilities  of  false  alarm 
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Figure  3  Probability  of  detection  for  the  incoherent  sum  of  10  pulses  from  a  Siverling  1  target  in 

K-distributed  clutter  with  v  =  0.1,  for  different  probabilities  of  false  alarm 
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3.3  Target  Fluctuation  Models 

Figure  4  shows  the  probability  of  detection  for  the  Swerling  target  models  in  noise  following 
incoherent  integration  of  16  pulses.  Figure  4  reproduces  the  top  graph  in  Figure  8.17  from 
reference  [1], 
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Figure  4  Probability  of  detection  for  Swerling  target  models  in  noise  following  incoherent  integration 

of  16  pulses,  with  a  probability  of  false  alarm  of  10' 4 


3.4  Convergence  of  the  Integral  over  the  Clutter  Power 

Gaussian  quadrature  with  n  points  is  exact  for  an  integrand  in  the  form  of  a  polynomial  of 
order  2n  [11,  Section  4.5].  The  integrand  here  is  not  a  simple  polynomial,  so  the  number  of 
quadrature  points  required  for  a  sufficiently  accurate  result  needs  to  be  established.  Figure  5 
shows  the  probability  of  false  alarm  for  a  single  pulse  return  from  K-distributed  clutter  with 
y  =  0.1,  for  different  numbers  of  quadrature  points,  along  with  the  analytic  solution.  Figure  6 
shows  the  relative  error  in  the  probability  of  false  alarm  in  Figure  5.  It  is  noticeable  that  the 
number  of  points  required  to  reduce  the  relative  error  to  a  satisfactory  level  is  greater  for  low 
thresholds.  Figure  7  shows  the  probability  of  detection  for  a  single  pulse  return  from  a 
Swerling  2  target  in  K-distributed  clutter  with  v  =  0.1  and  a  probability  of  false  alarm  of  10"“, 
calculated  with  different  numbers  of  quadrature  points.  Figure  8  shows  the  absolute  error  in 
the  probability  of  detection  in  Figure  7,  estimated  as  the  difference  between  the  Pd  for  the  set 
number  of  quadrature  points  and  that  using  150  points.  The  absolute  error  with  10  quadrature 
points  is  less  than  10"2  for  all  signal  to  clutter  ratios. 
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Figure  5  Probability  of  false  alarm  for  a  single  pulse  return  from  K-distributed  clutter  with  v  =  0.1 

for  different  numbers  of  quadrature  points.  The  infinite  limit  is  the  analytic  solution 
equation  (24). 


Figure  6  Relative  error  in  the  probability  of  false  alarm  in  Figure  5 
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Figure  7  Probability  of  detection  for  a  single  pulse  return  from  a  Swerling  2  target  in  K-distributed 

clutter  with  v  =  0.1  and  a  probability  of  false  alarm  of  10~2,  calculated  with  different 
numbers  of  quadrature  points 
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Figure  8  Absolute  error  in  the  probability  of  detection  in  Figure  7 
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3.5  Rayleigh  Limit 

Figure  9  shows  the  probability  of  detection  for  the  incoherent  sum  of  10  pulses  from  a 
Swerling  2  target  in  K-distributed  clutter,  for  different  values  of  the  clutter  shape  parameter  v. 
As  noted  in  Section  2.3,  the  largest  value  of  v  for  which  the  integration  over  the  clutter  power 
can  be  done  in  Matlab  is  171.  Figure  9  shows  this  is  not  quite  the  same  as  the  Rayleigh  limit  of 
v  — ►  go  where  the  clutter  can  be  treated  as  noise. 
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0.9 


0.8 


0.7 


0.6 
CL0  0.5 
0.4 
0.3 
0.2 
0.1 
0 

-5  0  5  10  15  20  25 

SCR  (dB) 

Figure  9  Probability  of  detection  for  the  incoherent  sum  of  10  pulses  from  a  Swerling  2  target  in 
K-distributed  clutter,  for  different  values  of  the  clutter  shape  parameter  v  and  a  probability 
of  false  alarm  ofl0~4 

3.6  Comparison  with  the  Method  of  Watts  and  Wicks 

Figure  10  shows  the  probability  of  detection  for  the  incoherent  sum  of  5  pulses  from  a 
Swerling  2  target  in  K-distributed  clutter,  for  different  values  of  the  clutter  shape  parameter  v. 
Figure  10  reproduces  Figure  6  from  reference  [3].  Calculation  of  182  data  points  in  Matlab 
required  0.11  s  using  the  method  of  Section  2  (the  curves  in  Figure  10)  and  0.04  s  using  the 
method  of  Watts  and  Wicks  [2,3]  (dots  in  Figure  10) 2.  The  speed  comparison  is  quite 
favourable  because  the  nodes  and  weights  for  the  numerical  integration  only  need  to  be 
calculated  once  for  each  curve  in  Figure  10.  In  a  simulation  model  the  clutter  shape  parameter 
will  normally  be  different  for  each  data  point,  so  the  nodes  and  weights  must  be  recalculated 
every  time.  In  this  situation  the  method  described  here  is  slower  than  the  method  of  Watts 


1 

10 

100 

170 

OO 


2  Calculations  carried  out  with  Matlab  R2009b  on  a  32  bit  desktop  computer  with  a  3.2  GHz  dual  core 
Pentium  processor,  running  Windows  XP. 
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and  Wicks  by  a  factor  of  about  30.  The  method  of  Watts  and  Wicks  is  an  approximation  which 
involves  only  linear  interpolation  in  a  table,  following  some  preliminary  calculations.  The 
method  described  here  is  slower,  but  more  accurate  and  more  flexible  in  that  it  allows  a  wider 
range  of  target  fluctuation  models  with  no  additional  work.  The  coefficients  required  for  the 
method  of  Watts  and  Wicks  have  only  been  published  for  Swerling  0,  1  and  2  targets  [3]; 
calculation  of  the  necessary  coefficients  for  other  target  models  would  be  quite  laborious. 
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Figure  10  Probability  of  detection  for  the  incoherent  sum  of  5  pulses  from  a  Swerling  2  target  in 
K-distributed  clutter,  for  different  values  of  the  clutter  shape  parameter  v  and  a  probability 
of  false  alarm  ofl0A.  The  curves  are  calculated  using  the  method  described  in  Section  2,  and 
the  dots  are  calculated  using  the  method  of  Watts  and  Wicks  [2,3]. 

3.7  Interpolation  Scheme 

Faster  evaluation  is  highly  desirable  for  simulation  models  used  in  operations  research.  The 
method  of  calculation  described  here  could  be  used  to  create  interpolation  tables  for  Pd,  which 
may  assist  in  speeding  up  computation  in  some  situations.  A  possible  interpolation  scheme  is 
explored  in  Appendix  C. 
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4.  Possible  Extensions 

4.1  Constant  False  Alarm  Rate  Processing 

Constant  false  alarm  rate  (CFAR)  processing  is  often  used  to  adapt  the  detection  threshold  to 
the  changing  clutter  background.  A  variety  of  schemes  are  used  [1,  Chapter  9],  some  more 
amenable  to  mathematical  analysis  than  others.  The  limiting  case  is  'ideal  CFAR'  in  which  the 
threshold  follows  the  clutter  background  exactly.  Ideal  CFAR  can  be  modelled  by  adapting 
the  calculations  described  above,  but  the  resulting  probability  of  detection  is  unrealistically 
high,  particularly  for  spiky  clutter  (small  v)  [1,  Section  9.2.3. 7].  The  fixed  threshold  calculation 
is  more  useful  for  operations  analysis  because  it  gives  an  estimate  of  worst  case  performance. 
Another  approach  is  to  calculate  the  detection  performance  for  a  fixed  threshold  and  then 
apply  a  separately  calculated  CFAR  loss  [1].  Shnidman  [9]  has  devised  mathematical 
expressions  for  cell  averaging  CFAR,  but  these  are  based  on  a  non-central  chi-square 
distribution  for  the  clutter  [13]  instead  of  the  K  distribution. 

4.2  KK  Distribution 

The  KK  distribution  is  a  linear  combination  of  two  K  distributions  with  different  parameters, 
used  when  non-Rayleigh  white  caps  or  sea  spikes  extend  the  tail  of  the  distribution  [14,15].  It 
should  be  straightforward  to  extend  the  method  of  calculation  described  here  to  the 
KK  distribution,  simply  by  replacing  equation  (1)  with  the  sum  of  two  integrals,  one  for  each 
component  of  the  distribution. 

4.3  Clutter  Correlations 

The  method  described  by  Ward,  Tough  and  Watts  [1]  allows  for  pulse  to  pulse  correlation  in 
the  clutter  speckle  with  an  effective  number  of  looks  L  which  is  less  than  N,  the  number  of 
pulses  integrated.  In  the  case  of  fixed  frequency  operation,  where  the  clutter  speckle  is 
constant  from  pulse  to  pulse,  L  =  1.  Military  radars  generally  use  frequency  agility  to 
decorrelate  the  clutter  from  pulse  to  pulse.  Residual  correlation  in  the  clutter  can  be  accounted 
for  with  the  effective  number  of  looks  L  <  N,  with  L  not  necessarily  an  integer  [16].  The  basic 
framework  of  the  calculation  is  similar,  but  the  clutter  speckle  is  included  with  the  target 
fluctuations  instead  of  adding  it  to  the  noise.  The  integration  over  the  local  clutter  power 
remains  the  same,  but  the  probability  of  detection  in  equations  (11)  and  (15)  is  replaced  with 
more  complicated  expressions  [1,  Chapter  8] 3. 

4.4  Target  Models 

Shnidman  [17]  has  further  expanded  on  the  Swerling  target  models  by  including  a  persistent 
component  in  the  target  reflection,  which  is  modelled  using  either  the  non-central  gamma  or 
the  non-central  gamma-gamma  distribution.  These  expanded  models  are  intended  for  targets 


3  The  author  has  coded  these  expressions  in  Matlab,  and  reproduced  Figures  8.7, 8.15  and  8.16  from  [1], 
but  the  code  is  rather  slow  and  inefficient.  It  is  not  clear  if  (or  how)  the  method  can  be  applied  in  the 
clutter  only  case,  when  the  noise  power  is  zero. 
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where  the  Swerling  cases  do  not  work  well  for  high  signal  to  noise  ratios,  and  the  Weibull  or 
log-normal  distributions  have  previously  been  used  to  model  the  signal  fluctuations.  In 
operations  analysis  we  are  mainly  concerned  with  finding  the  performance  limit  for  low 
signal  to  interference  ratios,  and  the  target  models  provided  by  the  gamma  distribution 
should  be  sufficient  for  this. 


5.  Conclusion 


A  method  for  calculation  of  the  probability  of  detection  in  K-distributed  clutter  and  noise  has 
been  described.  The  probability  of  detection  in  the  local  clutter  power  plus  noise  is  obtained 
from  a  series  solution,  and  this  solution  is  integrated  over  the  gamma  distribution  for  the  local 
clutter  power.  The  method  is  faster  than  direct  numerical  integration  over  the  K  distribution, 
and  it  may  be  suitable  for  use  in  simulation  models  for  operations  research.  It  is  slower,  but 
more  flexible  and  more  accurate,  than  the  method  of  Watts  and  Wicks,  an  approximate 
method  based  on  interpolation.  Faster  evaluation  is  highly  desirable  for  simulation  models 
used  in  operations  research.  The  method  described  here  could  be  used  to  create  tables  for 
interpolation,  which  should  produce  the  appropriate  mix  of  speed  and  accuracy  for 
operations  research  simulation  models. 
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Appendix  A:  The  Series  Expansion  for  the  Probability 

of  Detection  in  Noise 


Some  insight  into  the  nature  of  the  series  for  the  probability  of  detection  in  noise  can  be 
gained  using  the  Stirling  approximation  for  the  gamma  function  [18,  equation  (7.141)]: 


r(z)  =  (z  - 1) !  =  —  *  z  V  J— 

Z  V  Z 

or  lnT(z)  »  zlnz-z  +  |ln(2^)-|ln(z) 


(25) 


The  terms  in  the  first  sum  from  equations  (11)  and  (15)  are  in  the  form 


Hence 


Pm  =  ' 


Fm  -Y 

e 


Fm  -Y 

e 


ml  m f  (mi) 


In  pm  x  m In Y -  Y -  In m  - m In m  +  m-\ ln(2/r)  +  \ In (mi) 


=  m  In 


(Y\ 

\m) 


-Y  +  m-j  ln(2  Km) 


f 


•  m 


l-i 

1  2 
m 


Y 


2\ 


V 


— 1 

ym  j 


—  Y  +  m-\  ln(2^m) 


=  ~Mr-m)2  -jln(2xm) 


(26) 


(27) 


and 


exp(-2 


Km 


(28) 


so  for  large  Y  the  series  terms  form  a  normal  distribution  with  mean  Y  and  variance  m.  It 
follows  that  the  number  of  terms  required  to  obtain  a  desired  accuracy  in  Pd  is  proportional  to 
the  square  root  of  Y,  and  the  required  terms  are  centred  around  m  =  Y.  If  Y  is  large  the  terms 
are  very  small  for  small  m,  so  it  is  not  possible  to  start  the  series  evaluation  at  m  =  0,  because 
the  first  term  will  underflow  to  zero.  The  approximation  (28)  is  only  valid  for  large  Y  and  m. 
The  problem  of  finding  a  suitable  starting  value  for  m  is  addressed  by  Shnidman  [6]  using  the 
approximation 


pm  x  — J=exp(-«2)  where  u  =  J2Y - {m  - T)  -  ^ Jm  - \  (29) 

which  is  valid  for  wi  <  2 Y  and  m  »  Vi.  This  is  derived  from  an  approximation  for  the 
non-central  chi-square  distribution  [19]. 
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Appendix  B:  Chernoff  Bounds 

If  the  precision  of  the  computer  arithmetic  is  e,  Pd  evaluates  to  1  for  Pd  >  1  -  £.  Similarly,  Pd 
evaluates  to  zero  for  Pd  <  em,  where  em  is  the  smallest  number  which  can  be  represented  on  the 
computer.  For  the  double  precision  arithmetic  in  Matlab,  £  ~  10’15  and  em  ~  10  ~308.  Chernoff 
bounds  [9,  Section  IV  B],[20]  provide  a  means  to  avoid  unnecessary  calculations  in  these  cases. 

We  use  £  for  the  lower  bound  and  1  -  £  for  the  upper  bound. 

The  unit  step  function  can  be  bounded  from  above  by  exp(X(x-t)),  and  from  below  by 
1  -  exp (-A(x-f)),  for  positive  X  (Figure  11). 


Since 

/•OO  p  00 

J  f(x)dx  =  J  u(x  -  t)f(x)dx 


(30) 


we  can  bound  this  integral  from  above  or  below  by  replacing  the  unit  step  function  with  the 
appropriate  bound.  These  bounds  are  generally  easier  to  evaluate  than  the  original  integral, 
and  the  parameter  X  can  be  chosen  to  give  the  closest  bound. 


From  equations  (6)  and  (9),  the  probability  of  detection  for  N  radar  returns  from  a  fixed  target 
in  noise  is 
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The  Chernoff  bound  for  this  case  is 

CB{N,s,Y,A)  =  f  e 

J  0 


IS 

'  tL 

J  0 

U ) 

-(/l+s)  i 


exp (-AY  +  Nss 1/(1  -  A,)) 


(l-A)N 

The  value  of  A  which  gives  the  closest  bound  is  obtained  by  solving 


dHCB)  =  _y+ 


Ns  N 
-  + - =  0 


dA  (l-A)2  1  -A 

for  A.  Equation  (33)  can  be  rearranged  to  form  a  quadratic  in  A  with  the  solution 


A0=  1-  —  -. 
0  27 


v27y 


+  - 
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(32) 


(33) 


(34) 


where  the  negative  square  root  is  chosen  to  give  the  correct  sign  for  A o.  The  bounds  on  the 
probability  of  detection  are 


and 


PcJ[N,s,Y)<Cb[N,s,Y,A0 )  for  7  >s  +  N  ( A0  is  positive)  (35) 

Pd  ( N, s, 7)  >  1  - CB  ( N, s, 7,  A0 )  for  7  <s  +  N  (A0  is  negative).  (36) 


The  Chernoff  bound  Cb  is  calculated  for  each  set  of  parameter  values.  If  it  is  less  than  e,  the 
probability  of  detection  is  set  to  0  for  Y  >  s  +  N,  or  1  for  Y  <  s  +  N. 


The  probability  of  detection  for  a  fluctuating  target  in  noise  is 
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where  iFi(a;c;z)  is  the  confluent  hypergeometric  function  [4,  §9.210/1].  This  result  was 
obtained  by  Swerling  [8,  equation  (54)].  It  also  appears  as  equation  (9)  in  reference  [5],  but  the 
latter  version  contains  several  typographical  errors. 
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The  Chernoff  bound  for  this  case  is 
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Equation  (38)  is  evaluated  using  a  definite  integral  of  the  confluent  hypergeometric  function 
[4,  §7.621/5]: 


|  V"1  xFx{a\c\i)es,dt  =  T{c)s-c  (l-l/j)-0 


(39) 


The  value  of  A  which  gives  the  closest  bound  is  obtained  by  solving 
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for  A.  Equation  (40)  can  be  rearranged  to  form  a  quadratic  in  A  with  the  solution 

A0=b/2-jSj4- 
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Appendix  C:  Interpolation  Scheme 

In  a  calculation  of  radar  performance  as  a  function  of  ground  range,  the  signal  to  noise  ratio, 
the  clutter  to  noise  ratio  and  the  grazing  angle  all  vary  from  point  to  point.  The  clutter  shape 
parameter  v  also  changes,  because  it  depends  on  the  grazing  angle.  Thus  the  nodes  and 
weights  for  the  integration  over  the  clutter  power,  which  depend  on  v,  must  be  recalculated  at 
every  point.  This  makes  the  calculation  method  described  here  rather  slow  for  calculations  of 
radar  performance  such  as  those  in  [22],  It  may  be  preferable  to  tabulate  the  probability  of 
detection  as  a  function  of  signal  to  interference  ratio  (SIR) 4  for  different  values  of  v,  and  then 
use  interpolation  in  this  table  for  the  radar  performance  calculation. 

The  probability  of  detection  table  has  three  variable  parameters:  SIR,  CNR  and  v.  Variations  in 
many  parameters  of  interest,  such  as  sea  surface  wind  speed  and  direction,  and  platform 
altitude,  speed  and  heading,  affect  only  these  three  variables.  Other  radar  parameters  such  as 
the  false  alarm  number  and  the  number  of  pulses  integrated  remain  constant,  however  if  they 
are  changed  the  table  must  be  recalculated. 

Linear  interpolation  is  satisfactory  if  a  logarithmic  scale  is  used  for  the  variable  parameters, 
and  it  is  convenient  to  use  the  built-in  Matlab  function  for  interpolation  in  three  dimensions 
for  this  purpose.  Instead  of  interpolating  log(v),  it  may  be  advantageous  to  interpolate  the 
quantity 


77  =  l°gio(l  +  ^)  (43) 

The  SIR  for  Pd  =  0.5  is  approximately  linear  in  rj  (Figure  12),  and  the  Rayleigh  case  (v  — >  oo) 
corresponds  to  rj  =  0,  which  enables  interpolation  between  this  case  and  finite  values  of  v.  (The 
main  effect  of  varying  v  is  a  shift  of  the  Pd  curve  to  higher  signal  to  interference  ratios  as  v 
decreases,  which  is  quantified  by  the  SIR  required  for  Pd  =  0.5.)  Watts  and  Wicks  [2]  use  an 
empirical  expression  in  the  form  (43)  for  the  signal  to  clutter  ratio  for  Pd  =  0.5.  In  their 
expression  Vo  depends  on  the  false  alarm  number  and  the  target  fluctuation  model,  but  for  the 
purposes  of  interpolation  it  is  sufficient  to  use  a  constant.  If  we  take  Vo  =  9.9,  r\  ranges  from  0  to 
2  as  v  varies  from  infinity  down  to  0.1. 


4  SIR  =  SNR/  (1 +CNR)  in  clutter  and  noise;  SIR  =  SCR  in  clutter  only. 
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Figure  12  Signal  to  interference  ratio  for  Pd  =  0.5  plotted  against  the  parameter  rj  (equation  (43)). 

This  example  is  for  the  incoherent  sum  of  5  pulses  from  a  Swerling  2  target  in  K-distributed 
clutter,  with  a  probability  of  false  alarm  oflOA. 

The  effect  of  CNR  on  the  Pd  curve  is  highly  non-linear.  It  also  depends  on  v:  if  v  is  small, 
varying  CNR  has  a  significant  effect,  whereas  if  v  is  large  the  clutter  is  noise-like  so  varying 
CNR  has  little  effect  (Figure  13).  The  dependence  on  v  makes  it  difficult  to  linearise  the  effect 
of  CNR,  but  the  quantity 


#  =  k)g 


10 


1  +  ^— 

1  +  _J_ 

V  1  ^  CNR 


is  useful  for  interpolation.  As  CNR  varies  from  0  to  oo,  £  ranges  from  0  to  1. 


(44) 
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log1Q(1  +  9/(1  +  1/CNR)) 


Figure  13  Signal  to  interference  ratio  for  Pd  =  0.5  plotted  against  the  parameter  £  (equation  (44)),  for 
different  values  of  the  clutter  shape  parameter  v.  This  example  is  for  the  incoherent  sum  of  5 
pulses  from  a  Swerling  2  target  in  K-distributed  clutter,  with  a  probability  of  false  alarm  of 
10 

Steps  of  0.5  dB  in  SIR,  steps  of  0.125  in  tj,  and  steps  of  0.0625  in  £  were  found  to  be  sufficient  to 
create  a  table  for  successful  interpolation  of  the  probability  of  detection  for  fluctuating  targets. 
The  interpolation  table  contained  a  total  of  81x17x17  =  23409  points,  with  SIR  ranging  from  - 
5  dB  to  35  dB.  Many  points  in  the  interpolation  table  lie  outside  the  Chernoff  bounds 
(Appendix  B),  and  Pd  is  set  to  0  or  1  without  calculation  at  these  points.  The  nodes  and 
weights  for  integration  only  need  to  be  calculated  once  for  each  value  of  v  in  the  table,  and  the 
detection  threshold  only  needs  to  be  calculated  once  for  each  combination  of  v  and  CNR. 
Preparation  of  the  table  in  this  way  took  3.7  s,  compared  with  50.7  s  to  calculate  the  same  set 
of  Pd  values  individually,  with  the  nodes  and  weights  for  integration  and  the  detection 
threshold  computed  separately  for  every  point  5.  Clearly  interpolation  is  only  worthwhile  if 
creating  the  interpolation  table  takes  less  time  than  direct  calculation  of  Pd  at  the  required 
points.  The  Matlab  interpolation  functions  are  much  more  efficient  if  the  entire  matrix  of 
required  data  points  is  interpolated  with  a  single  function  call,  rather  than  interpolating  the 
points  one  by  one  in  a  nested  loop.  Interpolation  of  23120  points  took  0.01  s  with  a  single 
function  call. 


5  Calculations  carried  out  with  Matlab  R2010a  on  a  32  bit  desktop  computer  with  a  2.66  GHz  quad  core 
Pentium  processor,  running  Windows  XP. 
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The  Pa  curves  for  fixed  targets  can  be  very  steep,  so  a  finer  grid  is  needed  for  successful 
interpolation;  steps  of  0.1  dB  in  SIR,  steps  of  0.1  in  rj,  and  steps  of  0.0625  in  £  were  used,  for  a 
total  of  401x21x17  =  143157  points.  Computation  of  this  larger  table  took  12.3  s  (306  s  with  the 
nodes  and  weights  for  integration  and  the  detection  threshold  computed  separately  for  every 
value  of  Pd). 

The  number  of  variables  for  interpolation  can  be  reduced  to  two  if  v  is  replaced  by 

Veff  =  v(1+ ck)2  (45) 

In  this  case  the  calculation  is  done  for  clutter  only,  so  the  probability  of  detection  is  tabulated 
as  a  function  of  SCR  and  rj,  with  v  replaced  by  veff  in  equation  (43).  Equation  (45)  is  derived  in 
Appendix  D.  This  approximation  is  satisfactory  for  a  clutter  to  noise  ratio  of  ~10  dB  or  greater, 
or  for  any  CNR  if  v  is  -10  or  greater.  It  does  not  work  so  well  for  spiky  clutter  (small  v )  and 
low  CNR  [1,  Section  9.3. 3.1]. 
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Appendix  D:  Effective  Shape  Parameter 


The  effective  shape  parameter  veff  for  the  clutter  plus  noise  distribution  is  obtained  by 
matching  the  normalised  second  intensity  moment  of  the  combined  distribution  to  that  of  the 
K  distribution  for  clutter  alone  [21].  The  first  two  moments  of  the  clutter  plus  noise 
distribution  are  ([1],  pll3  with  b  =  v/pc): 

{Z)cn=Pc+Pn  (46) 


{z2  )cn  =  2  Pc  (  1  +  7)  +  4PcPn  +  2 P2, 


(47) 


For  the  clutter  only  distribution  we  have 


=  Pc 


=  2Pc{l  +  7) 


The  normalised  second  intensity  moment  is 


Hence 


(48) 

(49) 


(50) 


(51) 


For  the  clutter  plus  noise  distribution  we  have,  making  use  of  equations  (46)  and  (47), 


2  (z 


Veff  = 


=  _ 2{Pc+P„f _ 

2 Pc  (!  +  7)  +  4 Pc  P„  +  2P2n  -  2P)  -  4 P cPn  -  2Pn 
_  (1  +  ck) 

i+M 

=  ,/(1  +  ck) 


(52) 
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